source('../env.R')
It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2462 Columns: 10── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities_summary = communities %>% group_by(city_id) %>% summarise(
regional_pool_size = n(),
urban_pool_size = sum(relative_abundance_proxy > 0)
)
ggplot(communities %>% filter(relative_abundance_proxy > 0), aes(x = relative_abundance_proxy)) + geom_bar(stat = "bin")
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))
Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude data
community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv')) %>%
dplyr::select(city_id, mntd_normalised, fdiv_normalised, mass_fdiv_normalised, locomotory_trait_fdiv_normalised, trophic_trait_fdiv_normalised, gape_width_fdiv_normalised) %>%
left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
left_join(communities_summary) %>%
left_join(city_points[,c('city_id', 'city_nm')]) %>%
rename(city_name='city_nm') %>%
na.omit() %>%
arrange(city_id)
Rows: 341 Columns: 37── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (37): mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fdiv_normalised, fdiv_actual, fdiv_min, fdiv_max, fdiv_...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`
community_data_metrics
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))
Rows: 304 Columns: 5── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
fetch_normalised_traits = function(required_species_list) {
required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
required_traits$trophic_trait_normalised = normalise(required_traits$trophic_trait, min(required_traits$trophic_trait), max(required_traits$trophic_trait))
required_traits$locomotory_trait_normalised = normalise(required_traits$locomotory_trait, min(required_traits$locomotory_trait), max(required_traits$locomotory_trait))
required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), labels = c('Gape Width', 'Trophic Trait', 'Locomotory Trait', 'Mass'))
traits_normalised_long
}
fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))
Read in our phylogeny
phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')
Load resolve ecoregions
resolve = read_resolve()
Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees
to_species_matrix = function(filtered_communities) {
filtered_communities %>%
dplyr::select(city_id, jetz_species_name) %>%
distinct() %>%
mutate(present = TRUE) %>%
pivot_wider(
names_from = jetz_species_name,
values_from = "present",
values_fill = list(present = F)
) %>%
tibble::column_to_rownames(var='city_id')
}
community_nmds = function(filtered_communities) {
species_matrix = to_species_matrix(filtered_communities)
nmds <- metaMDS(species_matrix, k=2, trymax = 30)
nmds_result = data.frame(scores(nmds)$sites)
nmds_result$city_id = as.double(rownames(nmds_result))
rownames(nmds_result) = NULL
nmds_result
}
https://www.datacamp.com/tutorial/k-means-clustering-r
scree_plot = function(community_nmds_data) {
# Decide how many clusters to look at
n_clusters <- min(10, nrow(community_nmds_data) - 1)
# Initialize total within sum of squares error: wss
wss <- numeric(n_clusters)
set.seed(123)
# Look over 1 to n possible clusters
for (i in 1:n_clusters) {
# Fit the model: km.out
km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
# Save the within cluster sum of squares
wss[i] <- km.out$tot.withinss
}
# Produce a scree plot
wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
geom_point(size = 4) +
geom_line() +
geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
xlab('Number of clusters')
scree_plot
}
cluster_cities = function(city_nmds, cities_community_data, centers) {
set.seed(123)
kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
city_nmds$cluster = kmeans_clusters$cluster
cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}
plot_nmds_clusters = function(cluster_cities) {
cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}
plot_city_cluster = function(city_cluster_data_metrics, title) {
species_in_cluster = communities %>%
filter(city_id %in% city_cluster_data_metrics$city_id) %>%
dplyr::select(jetz_species_name, city_name, relative_abundance_proxy)
tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
gg_tree = ggtree(tree_cropped)
gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) +
geom_tile(aes(fill=relative_abundance_proxy)) +
scale_fill_gradientn(colours=c("#98FB98", "#FFFFE0", "yellow", "orange", "#FF4500", "red", "red"), values=c(0, 0.00000000001, 0.1, 0.25, 0.5, 0.75, 1), na.value = "transparent") +
theme_minimal() + xlab(NULL) + ylab(NULL) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(fill='Urban Proxy Abundance')
species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
gg_cities_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("FDiv") + ylim(0, 1)
gg_title = ggplot() + labs(title = title) + theme_minimal()
gg_presence %>% insert_top(gg_cities_mntd, height = 0.5) %>% insert_top(gg_cities_fd, height = 0.5) %>% insert_left(gg_tree, width = 0.75) %>% insert_right(gg_traits, width = 0.5) %>% insert_top(gg_title, height = 0.06)
}
REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive_using_abundance')
nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "San Jose" "Los Angeles" "Concord" "Tijuana" "Bakersfield"
[6] "Fresno" "Sacramento" "Mexicali" "Hermosillo" "Las Vegas"
[11] "Phoenix" "Tucson" "Durango" "Portland" "Chihuahua"
[16] "Aguascalientes" "Seattle" "Ciudad Juárez" "San Luis PotosÃ" "Mexico City"
[21] "Saltillo" "Vancouver" "Salt Lake City" "Albuquerque" "Monterrey"
[26] "Nuevo Laredo" "San Antonio" "Denver" "Austin" "Houston"
[31] "Dallas" "Oklahoma City" "Calgary" "New Orleans" "Kansas City"
[36] "Omaha" "St. Louis" "Bradenton" "Tampa" "Minneapolis [Saint Paul]"
[41] "Atlanta" "Orlando" "Louisville" "Chicago" "Indianapolis"
[46] "Milwaukee"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id))
Run 0 stress 0.1005678
Run 1 stress 0.1210166
Run 2 stress 0.1012776
Run 3 stress 0.1007176
... Procrustes: rmse 0.009402672 max resid 0.03496107
Run 4 stress 0.1000046
... New best solution
... Procrustes: rmse 0.007232701 max resid 0.0346979
Run 5 stress 0.1228731
Run 6 stress 0.1022505
Run 7 stress 0.1205421
Run 8 stress 0.1205421
Run 9 stress 0.1240951
Run 10 stress 0.1217438
Run 11 stress 0.1000046
... New best solution
... Procrustes: rmse 0.000004447461 max resid 0.00001302201
... Similar to previous best
Run 12 stress 0.1000046
... Procrustes: rmse 0.000003090941 max resid 0.00001098798
... Similar to previous best
Run 13 stress 0.121519
Run 14 stress 0.1012776
Run 15 stress 0.1234703
Run 16 stress 0.1217145
Run 17 stress 0.1007176
Run 18 stress 0.1217145
Run 19 stress 0.1205423
Run 20 stress 0.124095
*** Best solution repeated 2 times
nearctic_cities_nmds
scree_plot(nearctic_cities_nmds)
nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(nearctic_cities)
nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
ggplot() +
geom_sf(data = nearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neartic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neartic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neartic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neartic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
nearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neartic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster5.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Culiacán" "Guadalajara" "Morelia" "Acapulco"
[5] "Querétaro" "Cuernavaca" "Puebla" "Oaxaca"
[9] "Xalapa" "Veracruz" "Tuxtla Gutiérrez" "Villahermosa"
[13] "Guatemala City" "San Salvador" "San Pedro Sula" "Mérida"
[17] "Tegucigalpa" "Managua" "San José" "Cancún"
[21] "Guayaquil" "Chiclayo" "Panama City" "Trujillo"
[25] "Quito" "Havana" "Cali" "Lima"
[29] "Pereira" "Miami" "MedellÃn" "Ibagué"
[33] "Cartagena" "Kingston" "Bogota" "Barranquilla"
[37] "Bucaramanga" "Cúcuta" "Maracaibo" "Arequipa"
[41] "Barquisimeto" "Santo Domingo" "Maracay" "El Alto [La Paz]"
[45] "Caracas" "Cochabamba" "Viña del Mar [ValparaÃso]" "RÃo Piedras [San Juan]"
[49] "Barcelona" "Concepción" "Santiago" "Mendoza"
[53] "Salta" "Cordoba" "Asuncion" "Buenos Aires"
[57] "La Plata" "Ciudad del Este" "Montevideo" "Mar del Plata"
[61] "Porto Alegre" "São Paulo" "Santos" "Sao Jose dos Campos"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id))
Run 0 stress 0.134619
Run 1 stress 0.1414255
Run 2 stress 0.1348046
... Procrustes: rmse 0.0106512 max resid 0.05474793
Run 3 stress 0.1348237
... Procrustes: rmse 0.01071367 max resid 0.05473745
Run 4 stress 0.1405639
Run 5 stress 0.1414222
Run 6 stress 0.141222
Run 7 stress 0.1414222
Run 8 stress 0.1344509
... New best solution
... Procrustes: rmse 0.006925358 max resid 0.04568131
Run 9 stress 0.1402357
Run 10 stress 0.1348046
... Procrustes: rmse 0.006662375 max resid 0.04608377
Run 11 stress 0.1405635
Run 12 stress 0.1405642
Run 13 stress 0.134619
... Procrustes: rmse 0.006938373 max resid 0.04565142
Run 14 stress 0.1346368
... Procrustes: rmse 0.006829873 max resid 0.04569156
Run 15 stress 0.1405644
Run 16 stress 0.1348237
... Procrustes: rmse 0.006540716 max resid 0.04604158
Run 17 stress 0.134433
... New best solution
... Procrustes: rmse 0.001186138 max resid 0.006984211
... Similar to previous best
Run 18 stress 0.1405635
Run 19 stress 0.1344509
... Procrustes: rmse 0.001186071 max resid 0.006831159
... Similar to previous best
Run 20 stress 0.1344509
... Procrustes: rmse 0.001186504 max resid 0.006792118
... Similar to previous best
*** Best solution repeated 3 times
neotropic_cities_nmds
scree_plot(neotropic_cities_nmds)
neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(neotropic_cities)
neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
ggplot() +
geom_sf(data = neotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
neotropic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neotropic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Lisbon" "Porto" "Marrakesh" "Seville" "Dublin"
[6] "Málaga" "Madrid" "Glasgow" "Bilbao" "Liverpool"
[11] "Bristol" "Manchester" "Birmingham" "Leeds" "Newcastle upon Tyne"
[16] "Sheffield" "Nottingham" "Valencia" "London" "Toulouse"
[21] "Paris" "Barcelona" "Rotterdam [The Hague]" "Brussels" "Amsterdam"
[26] "Lyon" "Marseille" "Dusseldorf" "Nice" "Frankfurt am Main"
[31] "Zurich" "Oslo" "Stuttgart" "Hamburg" "Genoa"
[36] "Nuremberg" "Copenhagen" "Munich" "Berlin" "Dresden"
[41] "Rome" "Prague" "Stockholm" "Poznan" "Vienna"
[46] "Wroclaw" "Zagreb" "Gdansk" "Budapest" "Krakow"
[51] "Warsaw" "Helsinki" "Riga" "Belgrade" "Lviv"
[56] "Sofia" "Thessaloniki" "Saint Petersburg" "Minsk" "Athens"
[61] "Kyiv" "Istanbul" "Odesa" "Samsun" "Luxor"
[66] "Tel Aviv" "Jerusalem" "Tbilisi" "Yerevan" "Kuwait City"
[71] "Doha" "Abu Dhabi" "Dubai" "Bishkek"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id))
Run 0 stress 0.04961857
Run 1 stress 0.06079762
Run 2 stress 0.06489516
Run 3 stress 0.05546046
Run 4 stress 0.06348958
Run 5 stress 0.05395553
Run 6 stress 0.04967654
... Procrustes: rmse 0.02445083 max resid 0.1151427
Run 7 stress 0.06304717
Run 8 stress 0.07947601
Run 9 stress 0.05001033
... Procrustes: rmse 0.06742835 max resid 0.2224835
Run 10 stress 0.05392775
Run 11 stress 0.06304664
Run 12 stress 0.05258764
Run 13 stress 0.08298216
Run 14 stress 0.05660148
Run 15 stress 0.05105666
Run 16 stress 0.05672742
Run 17 stress 0.06304727
Run 18 stress 0.05479118
Run 19 stress 0.05101845
Run 20 stress 0.05471648
Run 21 stress 0.07604128
Run 22 stress 0.08017294
Run 23 stress 0.06304694
Run 24 stress 0.04961852
... New best solution
... Procrustes: rmse 0.0000829673 max resid 0.0003983299
... Similar to previous best
*** Best solution repeated 1 times
palearctic_cities_nmds
scree_plot(palearctic_cities_nmds)
palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 6)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(palearctic_cities)
palearctic_biomes = resolve[resolve$REALM == 'Palearctic',c('REALM')]
ggplot() +
geom_sf(data = palearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Palearctic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Palearctic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Palearctic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'))
Saving 14 x 6 in image
palearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Palearctic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Palearctic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'))
Saving 7.29 x 4.51 in image
palearctic_cities %>% filter(cluster == 6) %>% plot_city_cluster('Palearctic cluster 6')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Cape Town" "Johannesburg" "Pretoria" "Kigali" "Kampala" "Arusha" "Nairobi" "Addis Ababa" "Antananarivo"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id))
Run 0 stress 0.00009014786
Run 1 stress 0.00009691901
... Procrustes: rmse 0.0001926718 max resid 0.0004158943
... Similar to previous best
Run 2 stress 0.0003619111
... Procrustes: rmse 0.001981407 max resid 0.003127616
... Similar to previous best
Run 3 stress 0.00009161166
... Procrustes: rmse 0.0002339261 max resid 0.0003170711
... Similar to previous best
Run 4 stress 0.001912232
Run 5 stress 0.001215462
Run 6 stress 0.0006860128
Run 7 stress 0.00009535671
... Procrustes: rmse 0.00001612274 max resid 0.00002766371
... Similar to previous best
Run 8 stress 0.00009163152
... Procrustes: rmse 0.0002274177 max resid 0.0002943484
... Similar to previous best
Run 9 stress 0.001472905
Run 10 stress 0.3017338
Run 11 stress 0.00009434023
... Procrustes: rmse 0.0001876868 max resid 0.0004059187
... Similar to previous best
Run 12 stress 0.3017338
Run 13 stress 0.00009964034
... Procrustes: rmse 0.0004865478 max resid 0.001033375
... Similar to previous best
Run 14 stress 0.002511588
Run 15 stress 0.00009880586
... Procrustes: rmse 0.0001946382 max resid 0.0004211662
... Similar to previous best
Run 16 stress 0.00009893091
... Procrustes: rmse 0.0001949213 max resid 0.0004193287
... Similar to previous best
Run 17 stress 0.000468451
... Procrustes: rmse 0.002599256 max resid 0.004076805
... Similar to previous best
Run 18 stress 0.0009880508
Run 19 stress 0.00009746021
... Procrustes: rmse 0.0001852616 max resid 0.000383588
... Similar to previous best
Run 20 stress 0.00009997242
... Procrustes: rmse 0.000481814 max resid 0.001031423
... Similar to previous best
*** Best solution repeated 12 times
Warning: stress is (nearly) zero: you may have insufficient data
afrotropic_cities_nmds
scree_plot(afrotropic_cities_nmds)
afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 4)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(afrotropic_cities)
afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
ggplot() +
geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Afrotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Afrotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Afrotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster3.jpg'))
Saving 7.29 x 4.51 in image
afrotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Afrotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster4.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Srinagar" "Jamnagar" "Jammu" "Rajkot" "Bikaner" "Jodhpur"
[7] "Jalandhar" "Ahmedabad" "Bhavnagar" "Ludhiana" "Anand" "Udaipur"
[13] "Surat" "Vadodara" "Ajmer" "Chandigarh" "Vasai-Virar" "Mumbai"
[19] "Jaipur" "Delhi [New Delhi]" "Nashik" "Dehradun" "Kota" "Pune"
[25] "Haridwar" "Dhule" "Ujjain" "Indore" "Ahmadnagar" "Kolhapur"
[31] "Jalgaon" "Agra" "Aurangabad" "Sangli" "Belagavi" "Gwalior"
[37] "Budaun" "Bareilly" "Dharwad" "Bhopal" "Bhind" "Mangaluru"
[43] "Solapur" "Vijayapura" "Akola" "Latur" "Kannur" "Davanagere"
[49] "Thalassery" "Amravati" "Kalaburagi" "Kozhikode" "Guruvayur" "Malappuram"
[55] "Lucknow" "Thrissur" "Mysuru" "Kochi" "Alappuzha" "Nagpur"
[61] "Kollam" "Jabalpur" "Ettumanoor" "Hyderabad" "Coimbatore" "Bengaluru"
[67] "Thiruvananthapuram" "Tiruppur" "Faizabad" "Erode" "Prayagraj" "Pratapgarh"
[73] "Salem" "Dindigul" "Madurai" "Tiruchirappalli" "Durg" "Vellore"
[79] "Tirupati" "Raipur" "Bilaspur" "Vijayawada" "Puducherry" "Chennai"
[85] "Kathmandu" "Colombo" "Rajamahendravaram" "Patna" "Kandy" "Bihar Sharif"
[91] "Visakhapatnam" "Ranchi" "Brahmapur" "Jamshedpur" "Darjeeling" "Siliguri"
[97] "Cuttack" "Bhubaneshwar" "Jalpaiguri" "Berhampore" "Kolkata" "Krishnanagar"
[103] "Guwahati [Dispur]" "Agartala" "Silchar" "Dimapur" "Bangkok" "George Town"
[109] "Kuala Lumpur" "Phnom Penh" "Singapore" "Hong Kong" "Sha Tin" "Hsinchu"
[115] "Taichung" "New Taipei [Taipei]" "Tainan" "Denpasar" "Kaohsiung" "Kota Kinabalu"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id))
Run 0 stress 0.1190668
Run 1 stress 0.119963
Run 2 stress 0.1616726
Run 3 stress 0.1242016
Run 4 stress 0.1582327
Run 5 stress 0.1183841
... New best solution
... Procrustes: rmse 0.006657586 max resid 0.04580956
Run 6 stress 0.1219562
Run 7 stress 0.1153501
... New best solution
... Procrustes: rmse 0.02592979 max resid 0.2352578
Run 8 stress 0.1509741
Run 9 stress 0.155088
Run 10 stress 0.128717
Run 11 stress 0.1161599
Run 12 stress 0.1526359
Run 13 stress 0.1528153
Run 14 stress 0.1373095
Run 15 stress 0.1168773
Run 16 stress 0.152169
Run 17 stress 0.1296541
Run 18 stress 0.1256725
Run 19 stress 0.1526831
Run 20 stress 0.1364842
Run 21 stress 0.1153499
... New best solution
... Procrustes: rmse 0.00005138461 max resid 0.0003304756
... Similar to previous best
*** Best solution repeated 1 times
indomalayan_cities_nmds
scree_plot(indomalayan_cities_nmds)
indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 6)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(indomalayan_cities)
indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
ggplot() +
geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) +
geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 1) %>% plot_city_cluster('Indomalayan cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 2) %>% plot_city_cluster('Indomalayan cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 3) %>% plot_city_cluster('Indomalayan cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 4) %>% plot_city_cluster('Indomalayan cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 5) %>% plot_city_cluster('Indomalayan cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'))
Saving 7.29 x 4.51 in image
indomalayan_cities %>% filter(cluster == 6) %>% plot_city_cluster('Indomalayan cluster 6')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster6.jpg'))
Saving 7.29 x 4.51 in image